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Abstract 

The Birnbaum-Saunders regression model is commonly used in reliability studies. 
We address the issue of performing inference in this class of models when the number 
of observations is small. Our simulation results suggest that the likelihood ratio test 
tends to be liberal when the sample size is small. We obtain a correction factor which 
reduces the size distortion of the test. Also, we consider a parametric bootstrap 
scheme to obtain improved critical values and improved p-values for the likelihood 
ratio test. The numerical results show that the modified tests are more reliable in 
finite samples than the usual likelihood ratio test. We also present an empirical 
application. 

Key words: Bartlett correction; Birnbaum-Saunders distribution; Bootstrap; 
Likelihood ratio test; Maximum likelihood estimation. 



1 Introduction 



Different models have been proposed for lifetime data, such as those based 
on the gamma, lognormal and Weibull distributions. These models typically 
provide a satisfactory fit in the middle portion of the data, but oftentimes fail 
to deliver a good fit at the tails, where only a few observations are generally 
available. The family of distributions proposed by Birnbaum and Saunders 
(1969) can also be used to model lifetime data. It has the appealing feature 
of providing satisfactory tail fitting. This family of distributions was origi- 
nally obtained from a model in which failure follows from the development 
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and growth of a dominant crack. It was later derived by Desmond (1985) us- 
ing a biological model which followed from relaxing some of the assumptions 
originally made by Birnbaum and Saunders (1969). 

The random variable T is said to be Birnbaum-Saunders distributed with 
parameters a, r\ > 0, denoted B-S(a, 77), if its distribution function is given by 



where <£>(•) is the standard normal distribution function; a and 77 are shape 
and scale parameters, respectively. It is easy to show that r] is the median 
of the distribution: F T {rj) = $(0) = 1/2. For any k > 0, it follows that 
kT ~ B-S(a,krj). It is also noteworthy that the reciprocal property holds: 
T _1 ~ B-S(a, rj -1 ), which is in the same family of distributions [Saunders 



Rieck and Nedelman (1991) proposed a log-linear regression model based on 
the Birnbaum-Saunders distribution. They showed that if T ~ B-S(a,rj), 
then y = log(T) is sinh-normal distributed with shape, location and scale 
parameters given by a, fi = log(?7) and a = 2, respectively [y ~ SAf(a, //, a)]; 
see Section 2 for further details. Their model has been widely used and is an 
alternative to the usual gamma, lognormal and Weibull regression models; see 
Rieck and Nedelman (1991, § 7). Diagnostic tools for the Birnbaum-Saunders 
regression model were developed by Galea et al. (2004), Leiva et al. (2007) 
and Xie and Wei (2007), and Bayesian inference was developed by Tisionas 



Hypothesis testing inference is usually performed using the likelihood ratio 
test. It is well known, however, that the limiting null distribution (x 2 ) used 
in the test can be a poor approximation to the exact null distribution of the 
test statistic when the number of observations is small, thus yielding a size- 
distorted test; see, e.g., the simulation results in Rieck and Nedelman (1991, 
§ 5). Consider, for instance, the application in which interest lies in model- 
ing the die lifetime (T) in the process of metal extrusion, as in Lepadatu et 
al. (2005). As noted by the authors, the die life is mainly determined by its 
material properties and the stresses under load. They also note that the ex- 
trusion die is exposed to high temperatures, which can also be damaging. The 
covariates are the friction coefficient (xi), the angle of the die (x 2 ) and work 
temperature (2:3). Consider a regression model which also includes interaction 
effects, i.e., 

Vi = A) + Pixu + P2X21 + foxzi + f3iXiiX 2 i + fcxuxsi + (3 6 x 2 iX 3 i + e u 

where yi = log(Tj) and ~ SJ\f(a, 0, 2), i = l,2,...,n. There are only 15 
observations (n = 15), and we wish to test the significance of the interaction 




t > 0, 



(1) 



(1974)]. 



(2001). 
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effects, i.e., the interest lies in testing Ho '■ Pi — (3b — — 0. The likelihood 
ratio p- value equals 0.094, i.e., one rejects the null hypothesis at the 10% nom- 
inal level. Note, however, that the p-value is close to the significance level of 
the test and that the number of observations is small. Can the inference made 
using the likelihood ratio test be trusted? We shall return to this application 
in Section 6. 

The chief goal of our paper is to improve likelihood ratio inference in Birnbaum- 
Saunders regressions when the number of observations available to the prac- 
titioner is small. We do so by following two different approaches. First, we 
derive a Bartlett correction factor that can be applied to the likelihood ratio 
test statistic. The exact null distribution of the modified statistic is generally 
better approximated by the limiting null distribution used in the test than that 
of the unmodified test statistic. Second, we consider a parametric bootstrap 
resampling scheme to obtain improved critical values and improved p-values 
for the likelihood ratio test. 

The paper unfolds as follows. Section 2 introduces the Birnbaum-Saunders 
regression model. In Section 3, we derive a Bartlett correction to the likelihood 
ratio test statistic; we give a closed-form expression for the correction factor 
in matrix form. Special cases are considered in Section 4. Numerical evidence 
of the effectiveness of the finite sample correction we obtain is presented in 
Section 5; we also evaluate bootstrap-based inference. Section 6 addresses the 
empirical application introduced above (inferences on die lifetime in metal 
extrusion). Finally, concluding remarks are offered in Section 7. 



2 The Birnbaum Saunders regression model 

The density function of a Birnbaum-Saunders variate T is 



f T (t;a,rj) = 




where t, a, rj > 0. The density is right skewed, the skewness decreasing with a: 
see Lemonte et al. (2007, § 2). The mean and variance of T are, respectively, 



McCarter (1999) considered B-S{a,rj) parameter estimation under type II 
data censoring. Lemonte et al. (2007) derived the second order biases of the 
maximum likelihood estimators of a and 77, and obtained a corrected likelihood 




and 
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ratio statistic for testing hypotheses regarding a. Lemonte et al. (2008) pro- 
posed several bootstrap bias-corrected estimators of a and rj. Further details 
on the Birnbaum-Saunders distribution can be found in Johnson et al. (1995). 

The B-S(a,r]) survival function is Sxit) = 1 — Fxit), where Ft if) is given 
in (1). The hazard function is u(t) = /r(t)/S , 'r(t), where frit) is the corre- 
sponding density function. The hazard function u(t) equals zero at t — 0, 
increases up to a maximum value and then decreases towards a given posi- 
tive level; see Kundu et al. (2008). For a comparison between the Birnbaum- 
Saunders and lognormal hazard functions, see Nelson (1990). 

As noted in the previous section, Rieck and Nedelman (1991) showed that if 
T ~ B-S(a,rj), then y = log(T) follows a sinh-normal distribution with the 
following shape, location and scale parameters: a, \i = log(??) and a = 2, 
respectively, denoted y ~ SJ\f(a, /i, a). The density function of y is 

f(y;a,/i,a) = — -^=cosh^^— -^J expj-^sinh 2 j, y G JR. 

This distribution has a number of interesting and attractive properties [see 
Rieck (1989)]: (i) It is symmetric around the location parameter //; (ii) It is 
unimodal for a < 2 and bimodal for a > 2; (iii) The mean and variance of 
y are E(y) = /i and Var(y) = a 2 w(a), respectively. There is no closed-form 
expression for w(a), but Rieck (1989) obtained asymptotic approximations 
for both small and large values of a; (iv) If y a ~ SAf(a, fi,cr), then S a = 
2(y a — //)/ (acr) converges in distribution to the standard normal distribution 
when a — > 0. 

Rieck and Nedelman (1991) proposed the following regression model: 

yi = xJ(3 + Ei, i = l,2,...,n, (2) 

where yi is the logarithm of the ith. observed lifetime, xj = (xn,Xi2, ■ ■ ■ ,x ip ) 
contains the ith observation on the p covariates (p < n), (3 = (f3i, fa, ■ ■ ■ , {3 P ) T 
is a vector of unknown regression parameters, and £j ~ SAf(a, 0, 2). 

The log-likelihood function for a random sample y = (y±, . . . ,y n ) T from (2) 
can be written as 

n in 

£(6; y) = -- log(87r) + £ logfo) - - £ &, (3) 

i=l i=l 

where = (/3 T ,a) T , 

MO) = & = ^ cosh (S^)' ^ = ^ 2 = ^ sinh (S^) 



4 



and \X{ = xj(3, i — 1, 2, . . . , n. By differentiating (3) with respect to f3 r and a, 
we obtain 



de(o) i " r e. 



z2 




and 

ggg) n | 1 " 2 
(9a a %2 ' 

The score function for (3 can be written in matrix form as 

where X = (aq x 2 ■ ■ ■ x n ) T is the n x p design matrix (which is assumed to 
have full column rank) and s = s(0) is an n- vector whose ith element equals 

£il£i2 — 

Rieck and Nedelman (1991) obtained a closed-form expression for the maxi- 
mum likelihood estimator (MLE) of a 2 : 

where (3 is the MLE of (3. There is no closed-form expression for the MLE of 
(3. Hence, one has to use a nonlinear optimization method, such as Newton- 
Raphson or Fisher's scoring, to obtain (3. 1 

Let0 = T ,a) T be the MLE of = (f3 T ,a) T . Rieck and Nedelman (1991) 

showed that ~ Af p+1 (6, K(O)- 1 ), when n is large, ~ denoting approximately 
distributed; K(0) is Fisher's information matrix and K(6)~ l is its inverse. 
Also, K(0) is a block-diagnonal matrix given by K(9) = di&g{K((3), K a , a }: 
K{(3) = ipi(a)(X X)/4 is Fisher's information for (3 and K a>a = 2n/a 2 is the 
information relative to a. Also, 



M<*) = |l - erf j exp (J^ and ^i(«) = 2 + -i - ^^ ( a ), (4) 

erf (•) denoting the error function: 

2 



2 f x 2 

erf (x) = -= / e~* dt. 

\/7T JO 



1 All log-likelihood maximizations with respect to /3 and a in this paper were 
carried out using the BFGS quasi-Newton method with analytic first derivatives; 
see Press et al. (1992). The initial values in the iterative BFGS scheme were (3 = 
(X T X)~X T y for (3 and Va 2 for a, where a 2 is obtained from a 2 with/3 replaced 
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Details on erf (•) can be found in Gradshteyn and Ryzhik (2007). Since K(0) 
is block-diagonal, (3 and a are globally orthogonal [Cox and Reid (1987)] and 
f3 and a are asymptotically independent. It can be shown that when a is 
small, ipo(a) ~ a/\/2ir and ipi(ct) « 1 +4/a 2 ; when a is large, ipo(a) 1 and 
■0i(o;) ~ 2. 



3 An improved likelihood ratio test 



Consider a parametric model f(y;0) with corresponding log-likelihood func- 
tion £(6;y), where 6 = (0j,8 2 ) T is a Axvector of unknown parameters. The 
dimensions of B\ and 6 2 are k — q and g, respectively. Suppose the interest lies 
in testing the composite null hypothesis H - 2 = 2 ^ against H 2 : 2 ^ 6^ 2 \ 
where 6^ is a given vector of scalars. Hence, 0\ is a vector of nuisance pa- 
rameters. The log-likelihood ratio test statistic can be written as 

LR = 2{£(e;y)-£(0;y)}, (5) 

where = {O^Ojy and = (6? ,0 { 2 0)T ) T are the MLEs of = (0^,6>^) T 
obtained from the maximization of £(0; y) under 7ii and 7i , respectively. 

Bartlett (1937) computed the expected value of LR under 7i up to order n _1 : 
E(LR) = q + B(8) + 0(n~ 2 ), where B(0) is a constant of order 0{n' v ). It 
is possible to show that, under the null hypothesis, the mean of the modified 
test statistic 

TR LR 
LRb ~l + B(8)/q 

equals q when we neglect terms of order 0(n~ 2 ). The order of the approxima- 
tion remains unchanged when the unknown parameters in B(6) are replaced 
by their restricted MLEs. Additionally, whereas Pr(LR < z) — Pr(x^ < 
z) + 0(n~ r ), it follows that Pr(LR b < z) = Vi{x 2 q < z) + 0(n~ 2 ), a clear 
improvement. The correction factor c = 1 + B(0)/q is commonly refered to as 
the 'Bartlett correction factor'. 

Note that LR can be written as 

LR = 2{l(0; y) - £(6; y)} - 2{l{0 ] y) - £(6; y)}, 

where £(0; y) is the log-likelihood function at the true parameter values. Law- 
ley (1956) has shown that 



2E 



£(0;y)-£(0;y) = k + e k + 0(n ), 
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where e k is of order 0(n x ) and is given by 

^ ' {^rstu ^rstuvw) j (6) 

where Yl denotes summation over all components of 0, i.e., the indices r, s, t, u, v 
and w vary over all k parameters, and the A's are given by 

\ - K rs tuf K rstu (u) (su)\ 

A rstu — ro rv \ ^ ^rst ^ r°rt J ) 

A rstuvw = K TS K tU K VW ^K-rtv ^ g K stu^ (^) 

where « rs = E(d 2 £(0)/d6 r d6 s ), K rst = E(d 3 £(0)/d6 r d6 s d6 t ), 4} = dK rs /d6 t , 
etc., and — K rs is the (r,s) element of Fisher's information matrix inverse. 
Analogously, 



2E 



0;y)-£(0;y) = k-q + e k „ q + 0(n~ 2 ), 



where e k - q is of order 0(n~ r ) and is obtained from (6) when the sum J2' only 
covers the components of 0i, i.e., the sum ranges over the k — q nuisance 
parameters, since 6 2 is fixed under 7i - 

Under Tio, K(LR) = q + e k — t k ~ q + 0(n~ 2 ). Thus, it is possible to achieve 
a better Xq approximation by using the modified test statistic LR b = LR/c 
instead of LR, the Bartlett correction factor being c = 1 + B(0)/q, where 
B(0) = e k — e k - q . The corrected statistic LR b is \\ distributed up to order 
0{pT l ) under Ho- The improved test follows from the comparison of LR b and 
the critical value obtained as the appropriate x 2 q quantile. 

The corrected test statistic is usually written as LR b = LR/{1 + B(0)/q}. 
Nonetheless, there are alternative modified statistics that are equivalent to 
LR b to order O^ 1 ), such as LR* b = LRexp{-B(0)/q} and LR* b * = LR{1 - 
B(0)/q}. It is noteworthy that LR* h has an advantage over the other two 
specifications: it never assumes negative values. See Cribari-Neto and Cordeiro 
(1996) for further details on Bartlett corrections. 

In what follows, we shall derive the Bartlett correction factor for testing in- 
ference in the Birnbaum-Saunders regression model. The parameter vector 
is = ({3 T ,a) T , which is (p + l)-dimensional. Hence, we shall obtain e p+ i 
from (6), with the indices varying from 1 up to p + 1. 

Let Z = X(X T X)~ 1 X T = {zij} and Z d = diag{zn, z 22 , ■ ■ ■ , z nn }. Also, 
Z (2) = Z Z, zP = Z d Z d , etc., denoting the Hadamard (element- 
wise) product of matrices. We shall use the following notation for cumulants 
of log-likelihood derivatives with respect to (3 and a: U r = d£(0)/d[3 r , U a = 
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di{0)/da, U rs = d 2 £(6)/d(3 r d(3 s , U ra = d 2 £(0)/d[3 r da, U aa = d 2 £(0)/da 2 , 
U rst = d 3 £(0)/d[3 r d[3 s dp t , U rsa = d 3 £(0)/d[3 r d(3 s da, etc; n rs = E(U rs ), n ra = 
E(U ra ), K rst = E(U rst ), etc; = dK rs /d(3 t , k%°) = d 2 n ra /d(3 t da, etc. 

From the log-likelihood function in (3) we obtain the following cumulants: 



f O A / j — LI — 7 ~ I Lt "J — CtUt tj 

4 ^ « 

-n 2 + a 2 ^ lOn 

a J ~l cr 

3(2 + n 

5.1/; 

where 



K-rstu lf)2\0i) / % ir% is-E it-E iu i l^rsta 0, ^rsaa A / d Xj r Xj s , 

i=l a i=l 

^raaa and Kaaaa T 

or 



««) = -i{ 2 + ^-/l(^ + ^)*° (a) } 

and ^(oO an d i>i(cx) are defined in (4). For small a, we have ^(oO ~ —5/8 — 
1/a 2 ; for large a, ^(a) w —1/2. 

Using these cumulants and also making use of the orthogonality between f3 and 
a, we obtain, after long and tedious algebra (Appendix), e p+ i = e(a,p,X), 
where 

e(a,p,X)=e a (a,p)+ep(a,X), (8) 

with 

e a (a,p) = + 8 1 (a)p + 5 2 (a)p 2 ^j and 6/3(0;, X) = 5 3 (a)tr(zj 2) ). 

Here, tr(-) denotes the trace operator and 

r / x 2 + a 2 . . , a r / \ i 2 . . . 2a^ 3 (a) 1 

5o(a) = v^ao^' = 45o(a) \ 2T^ + 5o(a) - i^rr 

62(a) = 25 o (a0 2 , * s (a) = and ^(a) = ^ - ^ [1 + ±W a ). 

In expression (8) - our main result - we write e p+ i as the sum of two terms, 
namely e a (a,p) and e / g(a, X). The quantity €p(a, X) is obtained from (6) with 
X/ ranging over the components of (3, i.e. as if a were known. The quantity 
e a (a,p) is the contribution yielded by the fact that a is unknown (see the 
Appendix). Note that e a (a,p) depends on the design matrix only through its 
rank. More specifically, it is a second degree polynomial in p divided by n. 
Hence, e a (a,p) can be non-negligible if the dimension of f3 is not considerably 
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smaller than the sample size. It is also noteworthy that e(a,p,X) depends 

on a but not on (3. The dependency of e(a,p,X) on a occurs through 5i(a), 

S 2 (a) and 63(a). For small a, we have Si(a) ~ 1, 6 2 (a) ~ 1/2 and $3(0:) ~ 0. 

(2) 

For large a, 61(a) « 1, 62(a) « 1/2 and 63(a) —1/2. Furthermore, tr(Z<] ; ) 
establishes the dependency of e(a,p, X) on X. In other words, e(a,p, X) 
depends on the sum of squares of the diagonal elements of the hat matrix Z. 
In particular, if p — 1, i.e. if X has a single column, x = (x±, . . . ,x n ) T say, 

then ti(Z^) = J2?=i "A I (S?=i xf) , the sample kurtosis of x. 

Finally, it should be noted that expression (8) is quite simple and can be easily 
implemented into any mathematical or statistical/econometric programming 
environment, such as R [R Development Core Team (2006)], Ox [Cribari-Neto 
and Zarkos (2003); Doornik (2006)] and MAPLE [Abell and Braselton (1994)]. 



4 Special cases 

In this section we present closed-form expressions for the Bartlett correction 
factor in situations that are of particular interest to practitioners. The simpli- 
fied expressions are obtained from our more general result given in (8). 

At the outset, we consider the test of H : a = a^ against Hi- a 7^ a^°\ where 
is a given positive scalar and f3 is a vector of nuisance parameters. The 
Bartlett correction factor becomes c = 1 + B(6), where B(6) = e(a,p,X) — 
e@(a, X), and hence, B(6) = e a (a,p). Note that the correction factor depends 
on X only through its rank, p. In particular, when p = 1 (i.i.d. case), we have 



This formula corrects eq. (14) in Lemonte et al. (2007), which is in error. For 
small and large values of a, we have B(6) « ll/(6n). 

Oftentimes practitioners wish to test restrictions on a subset of the regression 
parameters. For instance, one may want to test whether a given group of co- 
variates are jointly significant. To that end, we partition (3 as (3 = ((3 1 ,/3 2 ) T , 
where (3 1 = (ft, (3 2 , . . . , P P - q ) T and (3 2 = (ft>-«+i> Pp-q+2, ■ ■ ■ , P P V are vec- 
tors of dimensions (p — q) x 1 and 5 x 1, respectively, and consider the test 
of Ho- f3 2 = (3^ against Hi : (3 2 7^ (3^ , where (3^ is a q- vector of known 
constants. The most common situation is that in which (3 2 ^ = 0. Note that 
(3 l and a are nuisance parameters. In accordance with the partition of (3, 
we partition X as X = (Xi X 2 ), where the dimensions of Xi and X 2 are 
n x (p — q) and n x q, respectively. The correction factor is c = 1 + B(0)/q, 




9 



where B(0) = e(a,p, X) — e(a,p — q, Xi). It is easy to obtain 

B{0) = -{ + S 2 (a)q(2p - q)\ + 5 3 (a)tr(Zf - zjj), 

with Zi = X^X^X^X^ = {z Uj } and Z lci = diag{z m , z 122 , z lnn }. 

Next, suppose we wish to test Hq. (3 = (3^ against Hi. (3 ^ (3^\ where (3^ 
is a p-vector of known constants and a is a nuisance parameter. The Bartlett 
correction factor is c = 1 + B(0)/p with B(9) = e(a,p, X) — e a (a, 0), which 
yields 

B(6) = -{5i(a)p + 5 2 (a)p 2 }+5 3 (a)ti(zP). 



5 Numerical evidence 

We shall now report Monte Carlo evidence on the finite sample performance of 
three tests in Birnbaum-Saunders regressions, namely: the likelihood ratio test 
(LR), the Bartlett-corrected likelihood ratio test (LR b ), and an asymptotically 
equivalent corrected test (LRl). 2 The model used in the numerical evaluation 
is 

Vi = (hxn + I3 2 x i2 H h PpX ip + e h 

where Xn = 1 and Si ~ SN{ot, 0, 2), % = l,2,...,n. The covariate values 
were selected as random draws from the U(0, 1) distribution. The number of 
Monte Carlo replications was 10,000, the nominal levels of the tests were 7 
= 10%, 5% and 1%, and all simulations were carried out using the Ox matrix 
programming language (Doornik, 2006). 

Table 1 presents the null rejection rates (entries are percentages) of the three 
tests. The null hypothesis is Ho : /3 p _i = f3 p = 0, which is tested against 
a two-sided alternative, the sample size is n = 30 and a = 0.5. Different 
values of p were considered. The values of the response were generated using 
Pi=P 2 = --- = P P - 2 = l- 

Note that the likelihood ratio test is considerably oversized (liberal), more so 
as the number of regressors increases. For instance, when p = 8 and 7 = 10%, 
its null rejection rate is 18.78%, i.e., nearly twice the nominal level of the test. 
The two corrected tests are much less size distorted. For example, their null 
rejection rates in the same situation were 11.82% (LR b ) and 11.13% (LRl). 

The results in Table 2 correspond to a — 0.5 and p — 6. We report results 
for samples sizes ranging from 20 to 200. The null hypothesis under test is 

2 We do not report results relative to LR^* since they were very similar to those 
obtained using LR\. 
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Table 1 



Nul 


. rejection rates; 


a = 0.5, n = 


30. 












7 = 10% 


7 = 5% 


7 = 1% 


P 


LR 


LR h 


LRT 


LR 


LRu 


LRT 


LR 


LR h 


LRT 


3 


12.69 


10.36 


10.22 


6.51 


4.98 


4.90 


1.75 


1.25 


1.23 


4 


13.44 


10.27 


10.07 


7.46 


5.41 


5.32 


1.90 


1.10 


1.04 


5 


14.77 


10.74 


10.45 


8.25 


5.53 


5.31 


2.21 


1.18 


1.14 


6 


15.94 


11.07 


10.53 


9.14 


5.55 


5.23 


2.54 


1.24 


1.17 


7 


17.28 


11.55 


10.88 


10.13 


5.69 


5.42 


2.95 


1.29 


1.19 


8 


18.78 


11.82 


11.13 


11.15 


6.44 


5.83 


3.38 


1.48 


1.31 


9 


19.92 


12.11 


11.00 


12.00 


6.33 


5.66 


3.82 


1.49 


1.25 



7^o : fa = f3% = 0. The figures in this table show that the null rejection rates of 
all tests approach the corresponding nominal levels as the sample size grows, 
as expected. It is also noteworthy that the likelihood ratio test displays liberal 
behavior even when n = 100. Overall, the corrected tests are less size distorted 
than the unmodified test. For example, when n = 50 and 7 = 5%, the null 
rejection rates are 7.49% (LR), 5.32% (LR b ) and 5.17% (LR* b ). 

Table 2 



Null rejection rates; a = 0.5, p = 6 and different sample sizes. 





7 = 10% 


7 


= 5% 




7 = 1% 


n 


LR 


LRb 


LRl 


LR 


LRb 


LRl 


LR 


LRb 


LRl 


20 


19.54 


12.04 


11.08 


11.97 


6.54 


5.87 


4.05 


1.58 


1.38 


30 


15.94 


11.07 


10.53 


9.14 


5.55 


5.23 


2.54 


1.24 


1.17 


40 


13.57 


10.14 


9.97 


7.45 


4.99 


4.81 


1.79 


1.03 


1.01 


50 


13.36 


10.72 


10.51 


7.49 


5.32 


5.17 


1.51 


1.02 


0.99 


100 


11.86 


10.46 


10.44 


5.90 


4.92 


4.88 


1.25 


1.04 


1.03 


200 


10.92 


10.14 


10.12 


5.57 


5.07 


5.07 


1.04 


0.96 


0.96 



Figure 1 plots relative quantile discrepancies against the associated asymptotic 
quantiles for the three test statistics. Relative quantile discrepancies are de- 
fined as the difference between exact (estimated by Monte Carlo) and asymp- 
totic quantiles divided by the latter. Again, p = 6 and we test the exclusion 
of the last two covariates. Also, n = 30 and a = 0.5. The closer to zero the 
relative quantile discrepancies, the more accurate the test. While Tables 1 and 
2 give rejection rates of the tests at fixed nominal levels, Figure 1 compares 
the whole distributions of the different statistics with the limiting null dis- 
tribution. We note that the relative quantile discrepancies of the likelihood 
ratio test statistic oscillates around 25% whereas for the two corrected statis- 
tics they are around 5% (LRh) and 3% (LRl). It is thus clear that the null 
distributions of the modified statistics are much better approximated by the 
limiting null distribution (x|) than that of the likelihood ratio statistic. 

Table 3 contains the nonnull rejection rates (powers) of the tests. Here, p = 
4, a = 0.5 and n = 30, 50, 100. Data generation was performed under the 
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0.00 - 



1 2 3 4 5 6 7 

asymptotic quantile 



Fig. 1. Relative quantile discrepancies plot: n = 30, p = 6 and a = 0.5. 

alternative hypothesis: (3 3 = (3 4 = 5, with different values of S (S > 0). We have 
only considered the two corrected tests since the likelihood ratio is considerably 
oversized, as noted earlier. Note that the two tests display similar powers. For 
instance, when n = 50, 7 = 5% and 5 = 0.5, the nonnull rejection rates are 
72.39% (LR b ) and 72.28% (LR* b ). We also note that the powers of the tests 
increase with n and also with 5, as expected. 

Table 4 presents the null rejection rates for inference on the scalar parameter 
a. Here, n = 30 and p = 2, 3 and 4. The null hypotheses under test are 
Ho- a = 0.5 and Ho: a = 1.0. The likelihood ratio test is again liberal. Note 
that the two corrected tests are much less size distorted. For instance, when 
p = 4, 7 = 5% and a = 1.0, the null rejection rates of the LR, LR^ and LRl 
tests were 12.03%, 5.20% and 4.02%, respectively. 

Our simulation results concerning tests on the regression parameters were 
obtained for a = 0.5. In practice, values of a between and 1 cover most 
of the applications; see, for instance, Rieck and Nedelman (1991). We shall 
now present simulation results for a wide range of values of a, namely a = 
0.1,0.3,0.5,0.7, 0.9,1.2,2,10,50 and 100. The new set of simulation results 
includes rejection rates of the likelihood ratio test that uses parametric boot- 
strap critical values (with 600 bootstrap replications). The parametric boot- 
strap can be briefly described as follows. We can use bootstrap resampling to 
estimate the null distribution of the statistic LR directly from the observed 
sample y = (yi, . . . , y n ) T . To that end, one generates, under Ho (i.e., im- 
posing the restrictions stated in the null hypothesis), B bootstrap samples 
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Table 3 

Nonnull rejection rates; a = 0,5, p = 4 and different samp le sizes. 



n 5 




LRl 


10% 5% 1% 


10% 5% 1% 


30 0.1 
0.2 
0.3 
0.4 
0.5 


13.20 6.91 1.57 
20.66 12.22 3.46 
33.07 21.63 7.49 
48 36 35 57 14 96 
65.11 51.59 26.42 


13.01 6.73 1.52 
20.40 12.02 3.30 
32.73 21.28 7.33 
48.08 35.20 14.61 
64.72 51.19 25.99 


50 0.1 
0.2 
0.3 
0.4 
0.5 


13.82 7.63 2.03 
25.89 16.03 5.11 
45.00 32.06 13.15 
65.07 52.09 28.18 
82.31 72.39 48.01 


13.71 7.60 1.99 
25.81 15.96 5.03 
44.86 31.95 13.07 
64.97 51.96 28.03 
82.14 72.28 47.88 


100 0.1 
0.2 
0.3 
0.4 
0.5 


18.66 11.02 2.91 
43.63 31.29 13.05 
73.47 62.39 37.40 
92.12 86.39 69.15 
98.76 97.34 89.98 


18.65 11.01 2.90 
43.61 31.28 13.02 
73.47 62.35 37.34 
92.11 86.37 69.11 
98.76 97.33 89.93 



Table 4 

Null rejection rates; inference on a; n = 30 and different values for p. 







Ho- 


a = 0.5 


Ho: 


a = 1.0 


p 




10% 


5% 


1% 


10% 


5% 


1% 


2 


LR 


12.76 


6.99 


1.82 


13.55 


7.26 


1.93 




LRb 


10.13 


5.15 


1.21 


10.52 


5.10 


1.08 




LRl 


9.90 


5.02 


1.20 


10.31 


4.94 


1.05 


3 


LR 


15.16 


8.64 


2.45 


16.10 


9.35 


2.70 




LRb 


10.46 


5.43 


1.16 


10.53 


5.06 


0.97 




LRl 


9.77 


5.02 


1.02 


9.65 


4.68 


0.81 


4 


LR 


18.16 


10.72 


3.39 


19.86 


12.03 


3.73 




LRb 


10.45 


5.37 


0.98 


10.73 


5.20 


0.75 




LRl 


9.29 


4.57 


0.72 


8.77 


4.02 


0.43 



(y* 1 , • • • , y* B ) from the assumed model with the parameters replaced by re- 
stricted estimates computed using the original sample (parametric bootstrap), 
and, for each pseudo-sample, one computes LR* b = 2{£(9* b ; y* b )-l(9* b ; y* b )}, 
b = 1,2, ...,B, where 9* b and 9* b are the maximum likelihood estimators 
of 9 obtained from the maximizations of £(9;y* b ) under Tl and Hi, re- 
spectively. The 1 — 7 percentile of LR* b is estimated by §i_ 7 , such that 
#{LR* b < qi^y/B = 1—7. One rejects the null hypothesis if LR > §!_ 7 . For 
a good discussion of bootstrap tests, see Efron and Tibshirani (1993, Chapter 
16). 

Figures in Table 5 provide important information. For all values of a the 
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Bartlett and the bootstrap corrections are very effective in pushing the rejec- 
tion rates toward the nominal levels. An advantage of the Bartlett correction 
over the bootstrap approach is that the first requires much less computational 
effort. It is noteworthy that as a grows, the rejection rates of the likelihood 
ratio test approaches the corresponding nominal levels, making the corrections 
less needed. 



Table 5 

Null rejection rates of Hp : §3 = §4 = 0; p = 4, n = 25 and different values for a. 







a = 


0.1 






a = 


0.3 




7 


LR 


LR b 


LRl 


LRboot 


LR 


LR b 


LRl 


LRboot 


10% 


14.33 


11.42 


11.32 


9.90 


14.45 


10.70 


10.44 


10.18 


5% 


7.88 


5.82 


5.74 


4.94 


8.01 


5.49 


5.22 


4.97 


1% 


1.95 


1.29 


1.23 


1.24 


2.07 


1.20 


1.13 


1.13 






a = 


0.5 






a = 


0.7 




7 


LR 


LR b 


LRl 


LRboot 


LR 


LR b 


LRl 


LRboot 


10% 


14.25 


10.23 


10.01 


10.23 


14.17 


10.61 


10.36 


10.14 


5% 


7.69 


5.17 


5.02 


5.12 


8.09 


5.35 


5.17 


5.28 


1% 


1.89 


0.91 


0.85 


1.24 


2.07 


1.12 


1.06 


1.02 






a = 


0.9 






a = 


1.2 




7 


LR 


LRb 


LRl 


LRboot 


LR 


LRb 


LRl 


LRboot 


10% 


13.96 


10.80 


10.60 


9.64 


13.49 


10.45 


10.29 


9.79 


5% 


8.03 


5.77 


5.61 


5.17 


7.51 


5.37 


5.26 


5.10 


1% 


2.29 


1.30 


1.26 


1.10 


1.90 


1.18 


1.16 


1.38 






a - 


= 2 






a = 


10 




7 


LR 


LRb 


LRl 


LRboot 


LR 


LR b 


LRl 


LRboot 


10% 


13.21 


10.87 


10.64 


10.01 


12.44 


11.23 


11.13 


9.75 


5% 


7.29 


5.61 


5.50 


5.08 


6.59 


5.81 


5.73 


4.80 


1% 


1.63 


1.08 


1.04 


1.21 


1.42 


1.17 


1.16 


0.98 






a = 


50 






a = 


100 




7 


LR 


LRb 


LRl 


LRboot 


LR 


LRb 


LRl 


LRboot 


10% 


11.26 


10.45 


10.43 


10.11 


10.87 


10.17 


10.12 


9.89 


5% 


5.60 


5.21 


5.20 


4.93 


5.67 


5.11 


5.07 


5.18 


1% 


1.19 


1.06 


1.06 


1.04 


1.23 


1.07 


1.07 


1.18 



We shall now try to shed some light on the issue of the possible effect of 
near-collinearity between the covariates X on the testing procedures. To do 
so, we performed an additional simulation experiment. We set p = 4 and 
selected the covariate values as follows: xn = 1, for % — 1, . . . , n, the values of 
x 2 were chosen as random draws from the U(0, 1) distribution and the pairs 
(373, 2^4) were selected as random draws from the bivariate normal distribution 
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A/*2(0, E), where the covariance matriz S has the following form 




The closer the value of p is to either extreme (—1 or 1), the stronger the 
linear relation between the covariates x 3 and x±. Table 6 presents simulation 
results for different values of p. The figures in this table suggest that the 
sample correlation between x 2 and x 3 does not have significant effect on the 
behaviour of the testing procedures. Hence, near-collinearity does not seem to 
a matter of concern. 

Table 6 

Null rejection rates of Tip : fa = @4 = 0; V = 4, n = 20 and different values for p. 







P = 


0.0 




7 


Li? 


LRb 


LR* b 


LRboot 


10% 


16.00 


11.15 


10.72 


10.20 


5% 


9.16 


5.33 


5.08 


4.90 


1% 


2.37 


1.27 


1.21 


1.16 






P = 


0.5 




7 


LR 


LR b 


lr; 


LRboot 


10% 


15.54 


10.72 


10.33 


10.14 


5% 


8.85 


5.73 


5.48 


5.22 


1% 


2.37 


1.15 


1.03 


1.10 






P = 


0.9 




7 


LR 


LRb 


LRl 


LRboot 


10% 


15.73 


11.14 


10.77 


10.31 


5% 


9.18 


6.06 


5.70 


5.40 


1% 


2.58 


1.15 


1.10 


1.21 



In all simulated situations, the likelihood ratio test was liberal. Of course, 
this is not a proof that this is always the case. Indeed, there may be situa- 
tions where it is conservative. Simulation results presented in the literature, 
however, suggest that the likelihood ratio test is often anti-conservative. For 
a theoretical justification in a simple situation, let random 
sample drawn from the N(p, a 2 ) distribution, with both p and a 2 unknown. 
Consider the test of Ho : p = po versus Tii : p ^ po- The asymptotic likeli- 
hood ratio test rejects 7i whenever LR > c 7 , where c 7 is the 1 — 7 quantile 
of the xi distribution. Equivalently, 7i is rejected when ^/n\~z — p \/a > 
fc (7,n), where z = YZ =1 Zj/n , a 2 = E™=i(^ _ z) 2 /{n - 1) and k(~/,n) = 
\J (exp(— c 7 /2) _2 / n — l)(n — 1). Table 7 shows the true levels of the likelihood 
ratio test, i.e. Pr(LR > c 7 ) evaluated at TIq, for different values of n and 7. 
Notice that, even in this simple situation, the likelihood ratio test is liberal 
when the sample is not large, in agreement with simulation results presented 
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elsewhere. See, for instance, Rieck and Nedelman (1991, Table 4) and Cordeiro 
et al. (1995). 

Table 7 

True level; normal distribution. 







7 




n 


1% 


5% 


10% 


5 


2.91 


9.79 


16.54 


8 


1.97 


7.64 


13.72 


12 


1.58 


6.64 


12.35 


20 


1.32 


5.93 


11.35 


50 


1.12 


5.36 


10.52 



6 An application 

We shall now turn to an empirical application that employs real data. We 
consider the investigation made by Lepadatu et al. (2005) on metal extrusion 
die lifetime. As noted by the authors (p. 38), "the estimation of tool life (fa- 
tigue life) in the extrusion operation is important for scheduling tool changing 
times, for adaptive process control and for tool cost evaluation." They also 
note (p. 39) that "die fatigue cracks are caused by the repeat application of 
loads which individually would be too small to cause failure." According to 
them, current research aims at describing the whole fatigue process by focus- 
ing on the analysis of crack propagation from very small initial defects. It is 
noteworthy that fatigue failure due to propagation of an initial crack was the 
main motivation for the Birnbaum-Saunders distribution. 

In Section 1, we explained that the interest lies in modeling the die lifetime 
(T) in the metal extrusion process, which is mainly determined by its material 
properties and by the stresses under load. The extrusion die is exposed to high 
temperatures, which can also be damaging. The covariates are the friction 
coefficient (xi), the angle of the die (x 2 ) and work temperature (£3). Consider 
a regression model which also includes interaction effects, i.e., 

Vi = A) + PlXli + fil^U + /3 3 X 3i + ^4X u X 2 i + P 5 XliX 3i + (3§X 2i X 3i + (9) 

where yi = log(Tj) and Si ~ SAf(a,0,2), i = 1,2, ... ,n. There are only 15 
observations (n = 15), and we want make inference on the significance of the 
interaction effects, i.e., we wish to test H ■ At = A = Pe = 0. The likelihood 
ratio test statistic (LR) equals 6.387 (p- value 0.094), and the two corrected 
test statistics are LR b = 4.724 (p- value 0.193) and LR* b = 4.492 (p- value 
0.213). The p- value of the bootstrap-based likelihood ratio test is 0.276. It is 
noteworthy that one rejects the null hypothesis at the 10% nominal level when 
the inference is based on the likelihood ratio test, but a different inference is 
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reached when the modified (Bartlett-corrected or bootstrap-based) tests are 
used. Recall from the previous section that the unmodified test is oversized 
when the sample is small (here, n = 15), which leads us to mistrust the 
inference delivered by the likelihood ratio test. 

We proceed by removing the interaction effects (as suggested by the three 
modified tests) from Model (9). We then estimate 

Hi = f3 + (3iXu + (3 2 X 2 i + P 3 X 3i + £i, 

i = 1, . . . , 15. The point estimates are (standard errors in parentheses): (3q = 
5.9011 (0.488), A = 0.7917(1.777), (3 2 = 0.0098(0.012), (3 3 = 0.0052(0.001) 
and a = 0.1982 (0.036). The null hypothesis Ho ■ P 3 = is strongly rejected 
by the four tests (unmodified and modified) at the usual significance levels. 
All tests also suggest the individual and joint exclusions of Xi and x 2 from the 
regression model. We thus end up with the reduced model 

Xji = Pq + p 3 x 3i + €i, 

i — 1, ... ,15. The point estimates are (standard errors in parentheses): P = 
6.2453 (0.326), % = 0.0052 (0.001) and a = 0.2039 (0.037). 

We now return to Model (9) and test Ho '■ Pi — Pi — Pi — Ph — P% — 
(exclusion of all variables but x 3 ). The null hypothesis is not rejected at the 
10% nominal level by all tests, but we note that the corrected tests yield 
considerably larger p-values. The test statistics are LR = 7.229, LR b = 5.610 
and LRl = 5.417, the corresponding p-values being 0.204, 0.346 and 0.367; the 
p-value obtained from the bootstrap-based likelihood ratio test equals 0.484. 



7 Conclusions 

We addressed the issue of performing testing inference in Birnbaum-Saunders 
regressions when the sample size is small. The likelihood ratio test can be 
considerably oversized (liberal), as evidenced by our numerical results. We 
derived modified test statistics whose null distributions are more accurately 
approximated by the limiting null distribution than that of the likelihood 
ratio test statistic. We have also considered a parametric bootstrap scheme to 
obtain improved critical values and accurate p-values for the likelihood ratio 
test. Our simulation results have convincingly shown that inference based on 
the modified test statistics can be much more accurate than that based on the 
unmodified statistic. The modified tests behave as reliably as a likelihood ratio 
test that relies on bootstrap-based critical values, with no need of computer 
intensive procedures. We recommend the use of the following statistics: LRj, 
and LRl. The latter has the advantage of only taking on positive values, which 
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is desirable. We have also presented an empirical application in which the use 
of the finite sample adjustment proposed in this paper can lead to inferences 
that are different from the ones reached based on first order asymptotics. 
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Appendix 



From (6), we have 

p+i p+1 

e p+l = ^ ^ -Vstu ^ ^ \rstuvw- 

r,s,t,u=l r,s,t,u,v,w=l 

Note that Ylr~s\ u=i ^rstu can be written as J2rs t «=i ^rstu plus terms in which at 
least one subscript equals a. It follows from the orthogonality between a and j3 that 
several terms equal zero. The non-zero terms are Ylr s=i ^rsaa, u=\ ^aatu and 
X aaaa . Also, J2 P riXu,v, w =i Rstuvw is given by YZ,,t,u,v,w=i Rstuvw plus the follow- 
ing terms: Er,s,t,tt=l ^rstuaa, Er,s,» ) w=l ^rsaavw, Et,u,w,M)=l ^aatuvw, Er,s=l ^rsaaaai 

Ylt^i^aatuaa, YX,w=i ^aaaavw and X aaaaaa . Here, we present the derivations of 
Ef s t u=i ^rstu and J2v w =i ^aaaavw The other terms can be obtained in a similar 
fashion. 

Note that Y%,s,t,u=i ^rstu = (l/^)J2r,s,t,u=i KrSKtUK rstu- Inserting the cumulants 
given in Section 3 into this expression we have 



p -\ v { n 

E ^ = 4 E V<2(a)E 

s,t,u=l r,s,t,u=l y i=l 



Xi r Xi s XitXi- 



i=l r,s,t,u=l 



P 



— ^E") E x irK rS Xi s ^ E 

i=i U,s=i J U,i*=i J 

= ^EK T ^)K T ^) = ^EK T ^) 2 > 

where xj = (xn,Xi2, ■ ■ ■ , Xi p ) represents the ith row of X and K& 3 = K(P)^ 1 = 
4(X T X)~ 1 /V'i(a) represents the inverse of Fisher's information matrix for (3. There- 
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fore, 



P A I I \ n 

^ rs ' M = ^R 2 2^1 ^ ( x x ) Xi t ■ 

r,s,t,u=l i=l 



Note that zu = xj (X T X) 1 Xi is the zth diagonal element of Zd given in Section 3. 
Hence, 

A _ 4V> 2 (q) A 2 _ 4 ^(a) f (2), 

r,stt=l M») 2 U Ma) 2 



From Yjv,w=l ^aaaavw = (l/4)(« aQ ) z 

K a«« Xyu,t«=i K K avw , we obtain 
A _ _cr^jm_ A vw \ 2 + a 2 ^ 1 

A «o«<™«< - 4n 2 2a 3 K ] a 3 X ™ X ™ f 

d,ui=1 I. i=l ) 



i=l \.u,mj=1 / i=l 



5(2 + a 2 ) ^ 5(2 + a 2 ) 5(2 + a 2 )p 



E J <+ I" « J , X 



2na 2 ipi(a) 2na 2 tp\(a) 2na 2 tp\(a) 
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